Damage and fluctuations induce loops in optimal transport networks 
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Leaf venation is a pervasive example of a complex biological network, endowing leaves with a 
transport system and mechanical resilience. Transport networks optimized for efficiency have been 
shown to be trees, i.e. loopless. However, dicotyledon leaf venation has a large number of closed 
loops, which are functional and able to transport fluid in the event of damage to any vein, including 
the primary veins. Inspired by leaf venation, we study two possible reasons for the existence of a 
high density of loops in transport networks: resilience to damage and fluctuations in load. In the 
first case, we seek the optimal transport network in the presence of random damage by averaging 
over damage to each link. In the second case, we seek the network that optimizes transport when 
the load is sparsely distributed: at any given time most sinks are closed. We find that both criteria 
lead to the presence of loops in the optimum state. 



PACS numbers: 89.75Hc, 89.75Da, 89.75Fb, 89.75.Kd 

Networks optimized for various transport characteris- 
tics are of evident importance in urban infrastructure 
and technology; more subtly, they are also of great im- 
portance in various natural settings, such as river net- 
works [Ij , or in biological settings where the network has 
been subjected to natural selection The question of 
convexity is central in the study of such networks; for 
many objective functions, the optimal network can be 
shown to be topologically a tree, a loopless graph such 
that cutting any bond disconnects it into two pieces. For 
example, in the case of a resistor network where the cost 
of a conductance C between two nodes is proportional 
to and the total cost of the network is bounded, the 
network which optimizes energy dissipation will be a uni- 
form "sheet" when 7 > 1, or a tree when 7 < 1, with a 
phase-transition-like change at 7 = 1 [3l [4] [18] . The net- 
works for 7 < 1 are related to optimal channel networks 
(OCNs) yij, widely used to study and model river net- 
works with large basins; the extensively studied case [5] 
is mostly devoted to arborescent distribution structures. 

Yet everyone is familiar with the feeling to the touch 
of decaying leaves on damp ground; they feel like a piece 
of gauze or cloth, exhibiting resistance when subject to 
shear forces. This would not be so were the skeleton of 
the leaves treelike, and in fact it evidently is not. Di- 
cotyledon leaf venation is one of many natural networks 
which contain loops, in fact recursively nested sets of 
loops [19]. Other examples include the vasculature of 
two-dimensional animal tissues such as the retina [6] |7] , 
the structural veins of some insect wings or the struc- 
tural weblike bracing of certain gorgonian corals such as 
sea fans, as well as the network of rivers at deltas or 
lowlands. Almost all artificial distribution networks con- 
tain loops, most evidently the streets in city plans (see 
Suppl. Materials). Animal vasculature of flat tissues, for 
instance retinas, shares many of the properties of leaf ve- 
nation (though mechanical stability is usually not one of 



its roles), while animal structural vein patterns, such as 
corals or insect wings, have mechanical roles unrelated 
to transport. As with any other natural system, we do 
not know a priori what (if any) is the functional being 
optimized in these examples, and our only starting point 
is to discard as a candidate any functional or objective 
function whose optimal networks differ dramatically in 
structure from the observed ones. Therefore we know 
from the outset that the functional analyzed in ^ could 
not possibly be the one being optimized in these cases, 
as for no value of 7 is the optimal network similar to the 
ones observed. 

Several studies have concentrated on the morphogen- 
esis of leaf venation patterns, trying to elucidate the 
mechanisms by which a dense set of loops is created 
[3 m [TOl [11]. In this work we take a complementary 
point of view, namely, we shall want to understand why 
such patterns were chosen by evolution in the first place. 
We make the ad hoc assumption that evolution of such 
patterns was not strongly constrained by choice of mor- 
phogenetic mechanism and that, were the patterns sub- 
optimal in some fashion, some other morphogenetic al- 
gorithm would have been evolved. The interplay be- 
tween morphogenetic mechanism and the optimization 
described herein is beyond the scope of this work. 

There are two important aspects of leaf physiology 
that bear directly on their survival ability. First, leaves 
are under constant attack, from the elements as well as 
pathogens, insects, and herbivores Were the leaf 

vascular network tree-like, damage to any vein would re- 
sult in the death of all the leaf section downstream from 
that vein, and hence entire leaf sectors would be observed 
to be dead. However, this is not the case, as illustrated 
in the experiment in Fig. [ij the loops in the venation 
permit flow to be routed around any injury to the veins, 
even in case of injury to the main vein. We thus seek to 
optimize leaf function in the presence of damage to the 



FIG. 1: Left: Flow routed around an injury in a lemon leaf. A circular cut (black circle) was made on the main vein of a 
lemon leaf; fluorescein (bright green/yellow) was injected post injury at the stem (lower left, below the field of view), and was 
observed to flow through the vein network around the injury, closing a number of loops and eventually reaching the tip of the 
leaf. The two secondary veins immediately "downstream" from the injury (i.e., the two attached immediately to the right of 
the black circle) are thus supporting reversed flow from the tertiary veins back to the primary vein. Full HD movie available 
as part of the Suppl. Materials. Right: Stomatal patchiness visualized by chlorophyll autofluorescence (455nm stimulation, 
broad red emission) in a lemon leaf. Entire domains of the leaf show all their stomata either closed (bright) or open (dark). 



veins. 

Second, the studies that found that a tree topology 
optimizes transport were carried out for the case of opti- 
mization under a spatially and temporally constant load. 
Yet in the case of, for instance, brain vasculature, we 
know that the system is sensitively optimized to deliver 
changing fluxes under constantly varying load; this is, 
in fact, the basis of the BOLD response that underlies 
functional magnetic resonance imaging. Indeed, there is 
a similar phenomenon in the case of leaves, called stom- 
atal patchiness. Stomatal patchiness refers to heteroge- 
neous stomatal aperture across the surface of the leaf 
blade (and thus water evaporation and photosynthetic 
activity [13]). Initially considered to be an artifact or 
epiphenomenon, stomatal patchiness is now believed to 
be part of plant physiology, a possible evolutionary adap- 
tation to situations in which ambient conditions would 
increase water flux beyond what it tolerable by the plant. 
A spatially and temporally irregular driving force in the 
form of stomatal patchiness would thus be contrary to 
the assumptions of the uniform load models described in 
Ref. [3l S] . Because of the nonlinearity of the optimiza- 
tion cost function, spatio-temporal irregularities can po- 
tentially alter the optima. In particular, the proofs of 
looplessness in Refs. O IH [14] all rely rather directly and 
essentially on constant fluxes. Therefore we also seek to 
optimize leaf function under conditions of varying load. 

Optimization under damage to veins implies, virtually 
by definition, the formation of loops. Even then, we 
shall find that more loops are formed than the minimum 
strictly required by topology, which is a ring joining all 



outermost nodes. More surprisingly, we shall find that 
optimization under a varying load leads to the formation 
of dense, recursively-looped structures. 

We now define the two models more precisely. We 
consider, following [3ii4|, a network consisting of nodes /c, 
joined by conductances Ckj] we shall denote by (j, k) the 
set of nodes j adjacent to a given node k. A conductance 
C is assumed to "cost" an amount C^, and the total cost 
of the network to be constant and equal to 1. 

At any given node k there are currents I^j through 
the edges of the graph, whose sum are net currents Xi 
which correspond either to evaporation through stom- 
ata or water injection at the stem. The currents are 
driven by differences in the potential Vi at the nodes, as 
hj = Ckj{Vk - Vj): 

^C,j{V,-Vj)=Ik (2) 

an equation whose inversion yields the Vk and conse- 
quently the Ikj. Finally, the functional being optimized 
is the total power dissipation 

p = lT.T.^M-yj)'' (3) 

These equations need to be supplemented with the stip- 
ulation that Xi = N — 1 (where k = 1 is the stem) and 
Xk = —1 for k > 1. 
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FIG. 2: Loops as a result of optimizing under damage to links 
(left column) and under a fluctuating load (right column) . In 
all plots the vein thickness is proportional to (7*^^+^/^^/^. 



In the first model (robustness to damage) , we compute 
the power dissipated when one of the links (e.g. ab) is 
broken: C^j = Ckj{l - Sakhj - ^ajhk)- For each set 
of conductivities with a cut bond we can compute V^^ 
and subsequently obtain the dissipation P^^ . The total 
dissipation is defined as: 



R 



Y^pa 

{ah) 



(4) 



which is the functional being minimized. Note that if 
breaking {ah) disconnects the graph, then P^^ becomes 
infinite. Thus, finiteness of R requires that breaking any 
one bond should not disconnect the graph. This can be 
satisfied by adding a perimeter ring of conductances (no 




FIG. 3: Normalized network dissipation vs 7. Red line: ho- 
mogeneous load, blue line: cut bond, black line: fluctuating 
sink. The number of network nodes was N — 949. Inset: 
Close up of the crossover for the fluctuating sink case. Solid 
line: medium size network [N — 949), dashed line: small 
network {N = 285), dotted line: big network {N = 1589). 



matter how small) to the terminal nodes of any tree. If 
the tree had T terminal nodes, then this would result in 
T facets in the graph. Many leaves have such a rim of 
veins [12]. However this only guarantees finiteness; to 
actually minimize the value of R more loops than that 
are required, as shown in the right column of Fig. [2] 

In our second model, we consider fluctuating load by 
introducing a single moving sink. We define = dok — 
5aki i-e., a single source at the stem and a single sink at 
node a. This induces new potentials ^ power loss 
and similarly to the cut bond model: 



(5) 



Optima of this functional are shown in Fig. [2j right col- 
umn. 

The model in [3l IH [14] was already computationally 
costly enough that it could only be optimized through 
relaxation techniques, or Monte Carlo methods that in- 
volved searching the space of tree graphs. The exten- 
sions presented in this paper cannot be optimized by 
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such a Monte Carlo method as the optima are not trees. 
We alternatively use an annealing method, described in 
the Supplementary Material. However, much like any 
method that does not use exhaustive search, the min- 
ima found cannot be guaranteed to be global minima. 
The method required the evaluation of the inverse of an 
N X N matrix G that depends on the conductivity. For 
computational efficiency in the cut bond model, rather 
than inverting a different matrix for each of the terms 
in the ab sum, we perturbatively expand G and show 
that it can be evaluated using only the inverse of the 
original matrix of the full network and a sparse matrix 
multiplication, dramatically speeding up the evaluation. 
To simplify the calculation in the fluctuating load model 
we similarly use the sparse nature of the X^. Details are 
given in Supplementary Material. 

Similarly to the results found in [3l [H [14] the cut bond 
model exhibits a transition from the 7 < 1 hierarchical 
structure case, to the 7 > 1 case, where the network is 
a uniform sheet. However, unlike the simple tree model, 
the veins anastomose and we observe the formation of 
nested loops, reminiscent of the ones seen in real leaves. 
Indeed, evolutionary trends in vascular plants indicate a 
tendency of the vascular system to develop redundancy 
and hierarchical network patterns [12 . 

As shown in Fig. [3j where we plot the ratio of the en- 
ergy of the network to the energy of a network with a 
constant conductivity distribution and same 7 (normal- 
ized dissipation), the results for the bond case are very 
similar to the uniform load case. They both exhibit a 
cusp at 7 = 1. The situation is drastically different for 
the moving sink model, where the phase transition has 
been replaced by a crossover shown in the inset of Fig. |3] 

Although in the sink model we still observe hierarchi- 
cal loops, the quantitative aspects of the the networks in 
the two cases are very different. For example, the average 
node valency is smaller for the sink model for 7 < 0.4, 
and becomes larger after that. In both cases, for 7 > 1 
it asymptotically approaches 6. Moreover, whereas for 
the broken bond model fractal-like features of the net- 
work persist until 7 < 1, for the sink model fractality 
disappears for 7 > 0.5 [15]. 

To conclude, it is widely yet incorrectly stated that 
many natural distribution networks, such as animal vas- 
culature or tree leaves, are treelike [16[|T7]; even the most 
cursory of visual inspections rapidly refutes such state- 
ments. Two fundamental misconceptions underlie these 
statements. First, an unstated yet universal implication 
is that the trees are hierarchically ordered, i.e., as they 
are inevitably depicted as regular binary trees, one gets 
the impression of an ordered branching from the top level, 
first order to second and subsequent orders. Second, the 
venous and arterial tree are described as touching each 
other strictly at the level of the "terminal nodes" of each 
tree, the capillaries, by implication the only member of 
the hierarchy where it is ambiguous whether it's a vein 



or an artery. This description is disproved by looking at 
real leaves or images of, for instance, retina, where it is 
evident that venous capillaries of all orders impinge di- 
rectly on primary arteries and viceversa, and that, more- 
over, veins of any order branch directly from higher-order 
veins of much larger order [6l [7] . This topological disor- 
der permits circulation in case of obstacles, but may also, 
as our results above suggest, confer superior deliverance 
to fluctuating loads. 

We're deeply indebted to Raymond Raad, Steffen 
Bohn, and Yves Couder; we'd like to acknowledge help 
from Keith Mott. During preparation of this manuscript 
we became aware of the work of Francis Corson. 
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The parameter 7 is related to the physical nature of the 
network in question; i.e. assuming equal length of links, 
for a pipe network with Poiseuille flow 7 = 0.25 whereas 
for a resistor network (uniform resistors, resistance only 
depending on thickness of wire) 7=1. 
Monocotyledon leaf venation also contains loops, al- 
though less evident. The prominent parallel veins are reg- 
ularly connected with smaller perpendicular veins. 



